Se pretende estimar \(\mathbb{E}_{f}[g(X_{1}, X_{2})]\), donde \[ f(x_{1}, x_{2}) \propto \begin{cases} \mathrm{e}^{-4 (x_{2} - x_{1}^{2})^{2} + (x_{2} - 1)^{2}} &\text{si } x_{2} \leq 2, \\ 0 &\text{en otro caso.} \end{cases} \qquad g(x_{1}, x_2) = \sqrt{x_{1}^{2} + x_{2}^{2}} \]
library(plotly)
f_u <- function(x1, x2) {
exp(-4 * (x2 - x1^2)^2 + (x2 - 1)^2) * (x2 <= 2)
}
malla <- data.frame(
x1 = (-250:250) / 100,
x2 = (-250:250) / 100
)
malla |>
plot_ly(x = ~x2, y = ~x1) |>
add_surface(z = outer(malla$x1, malla$x2, f_u))
Vamos a realizar la estimación requerida a partir de un paseo
aleatorio de Metropolis, haciendo uso del paquete mcmc.
Para ello, en primer lugar debemos definir el logaritmo de la densidad
objetivo como una función que toma como argumento un vector con los dos
valores \(x_{1}\) y \(x_{2}\).
log_f_u <- function(x) {
x1 <- x[[1]]
x2 <- x[[2]]
if (x2 <= 2) {
-4 * (x2 - x1^2)^2 + (x2 - 1)^2
} else {
-Inf
}
}
A continuación, iniciamos el paseo aleatorio de Metropolis a partir
del estado \(x_{1} = x_{2} = 0\), que
según el grÔfico se encuentra cerca de donde la densidad toma su mÔximo.
Podemos comprobar que la función metrop devuelve la
realización de la cadena de Markov en forma de matriz con dos columnas,
ya que cada estado es una tupla de dos nĆŗmeros reales.
library(mcmc)
paseo_Metropolis <- metrop(log_f_u, initial = c(0, 0), nbatch = 1e4)
head(paseo_Metropolis$batch)
## [,1] [,2]
## [1,] 0.0000000 0.0000000
## [2,] 0.0000000 0.0000000
## [3,] 0.0000000 0.0000000
## [4,] 0.7588121 -0.2884623
## [5,] 0.5888510 -0.2832625
## [6,] 0.5888510 -0.2832625
La regla empĆrica establece que para acelerar la convergencia de la cadena de Markov a la densidad objetivo, el porcentaje de estados aceptados deberĆa ser del 25 %.
paseo_Metropolis$accept
## [1] 0.2293
paseo_Metropolis <- metrop(paseo_Metropolis, scale = .9)
paseo_Metropolis$accept
## [1] 0.2662
paseo_Metropolis <- metrop(paseo_Metropolis, scale = .95)
paseo_Metropolis$accept
## [1] 0.2416
paseo_Metropolis <- metrop(paseo_Metropolis, scale = .93)
paseo_Metropolis$accept
## [1] 0.2471
Como diagnóstico de convergencia, los siguientes grÔficos parecen indicar que la cadena de Markov ha convergido a la densidad deseada.
library(coda)
paseo_Metropolis$batch |>
mcmc() |>
traceplot()
library(ggplot2)
ggplot(
as.data.frame(paseo_Metropolis$batch),
aes(x = V1, y = V2)
) +
geom_contour(
data = expand.grid(x1 = (-250:250) / 100, x2 = (-250:250) / 100),
mapping = aes(x = x1, y = x2, z = f_u(x1, x2))
) +
geom_point(size = .5)
El grĆ”fico de autocorrelación parece indicar que, al aplicar el mĆ©todo de las medias por lotes, la longitud de estos deberĆa ser de al menos 50 estados.
g <- function(x) {
x1 <- x[[1]]
x2 <- x[[2]]
sqrt(x1^2 + x2^2)
}
library(purrr)
paseo_Metropolis$batch |>
array_branch(margin = 1) |>
map_dbl(g) |>
mcmc() |>
autocorr.plot(lag.max = 100)
Estamos ya en condiciones de obtener la estimación requerida, asà como un intervalo de confianza para la misma.
paseo_Metropolis <-
metrop(paseo_Metropolis, nbatch = 1e3, blen = 1e2, outfun = g)
paseo_Metropolis$batch |>
mcmc() |>
autocorr.plot()
estimacion <- paseo_Metropolis$batch |>
array_branch(margin = 2) |>
map_dbl(mean)
error_estandar <- paseo_Metropolis$batch |>
array_branch(margin = 2) |>
map_dbl(var) |>
magrittr::divide_by(paseo_Metropolis$nbatch) |>
sqrt()
probabilidad_cobertura <- 0.95
alfa <- 1 - probabilidad_cobertura
percentil <- qnorm(1 - alfa / 2)
intervalo_confianza <- estimacion + c(-1, 1) * error_estandar * percentil
Se obtiene entonces una estimación de 0.8557277, con (0.8365232, 0.8749322) un intervalo de confianza con probabilidad de cobertura 0.95.